//==================================================== file = genzipf.c =====
//=  Program to generate Zipf (power law) distributed random variables      =
//===========================================================================
//=  Notes: 1) Writes to a user specified output file                       =
//=         2) Generates user specified number of values                    =
//=         3) Run times is same as an empirical distribution generator     =
//=         4) Implements p(i) = C/i^alpha for i = 1 to N where C is the    =
//=            normalization constant (i.e., sum of p(i) = 1).              =
//=-------------------------------------------------------------------------=
//= Example user input:                                                     =
//=                                                                         =
//=   ---------------------------------------- genzipf.c -----              =
//=   -     Program to generate Zipf random variables        -              =
//=   --------------------------------------------------------              =
//=   Output file name ===================================> output.dat      =
//=   Random number seed =================================> 1               =
//=   Alpha vlaue ========================================> 1.0             =
//=   N value ============================================> 1000            =
//=   Number of values to generate =======================> 5               =
//=   --------------------------------------------------------              =
//=   -  Generating samples to file                          -              =
//=   --------------------------------------------------------              =
//=   --------------------------------------------------------              =
//=   -  Done!                                                              =
//=   --------------------------------------------------------              =
//=-------------------------------------------------------------------------=
//= Example output file ("output.dat" for above):                           =
//=                                                                         =
//=   1                                                                     =
//=   1                                                                     =
//=   161                                                                   =
//=   17                                                                    =
//=   30                                                                    =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 genzipf.c                                                 =
//=-------------------------------------------------------------------------=
//=  Execute: genzipf                                                       =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (11/16/03) - Genesis (from genexp.c)                      =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <assert.h>             // Needed for assert() macro
#include <stdio.h>              // Needed for printf()
#include <stdlib.h>             // Needed for exit() and ato*()
#include <math.h>               // Needed for pow()
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/mman.h>

//----- Constants -----------------------------------------------------------
#define  FALSE          0       // Boolean false
#define  TRUE           1       // Boolean true
//----- Function prototypes -------------------------------------------------
int zipf(double alpha, int n);  // Returns a Zipf random variable
double rand_val(int seed);         // Jain's RNG

using namespace std;

//===== Main program ========================================================
int main(void) {
	//FILE   *fp;                   // File pointer to output file
	int fd, result;
	char file_name[256];        // Output file name string
	char temp_string[256];      // Temporary string variable
	double alpha;                 // Alpha parameter
	double n;                     // N parameter
	int num_values;            // Number of values
	int zipf_rv;               // Zipf random variable
	int i;                     // Loop counter
	int norm;					// Normalized
	char * map;
	int * iptr;
	float * fptr;

	// Output banner
	printf("---------------------------------------- genzipf.c ----- \n");
	printf("-     Program to generate Zipf random variables        - \n");
	printf("-------------------------------------------------------- \n");

	// Prompt for output filename and then create/open the file

	printf("Output file name ===================================> ");
	scanf("%s", file_name);
	fd = open(file_name, O_RDWR | O_CREAT | O_TRUNC, (mode_t) 0600);
	if (fd == -1) {
		printf("ERROR in creating output file (%s) \n", file_name);
		exit(1);
	}

	// Prompt for random number seed and then use it
	printf("Random number seed (greater than 0) ================> ");
	scanf("%s", temp_string);
	rand_val((int) atoi(temp_string));

	// Prompt for alpha value
	printf("Alpha value ========================================> ");
	scanf("%s", temp_string);
	alpha = atof(temp_string);

	// Prompt for N value
	printf("N value ============================================> ");
	scanf("%s", temp_string);
	n = atoi(temp_string);

	// Prompt for number of values to generate
	printf("Number of values to generate =======================> ");
	scanf("%s", temp_string);
	num_values = atoi(temp_string);

	// Prompt for normalization
	printf("Normalize output to values between 0 & 1 ==============> ");
	scanf("%s", temp_string);
	norm = atoi(temp_string);

	// Output "generating" message
	printf("-------------------------------------------------------- \n");
	printf("-  Generating samples to file                          - \n");
	printf("-------------------------------------------------------- \n");

	size_t FILESIZE = (num_values * (sizeof(int) + sizeof(float)));
	result = lseek(fd, FILESIZE - 1, SEEK_SET);
	if (result == -1) {
		close(fd);
		perror("Error calling lseek() to 'stretch' the file");
		exit(EXIT_FAILURE);
	}

	result = write(fd, "", 1);
	if (result != 1) {
		close(fd);
		perror("Error writing last byte of the file");
		exit(EXIT_FAILURE);
	}

	map = (char*) mmap(0, FILESIZE, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
	if (map == MAP_FAILED) {
		close(fd);
		perror("Error mmapping the file");
		exit(EXIT_FAILURE);
	}

	// Generate and output zipf random variables
	iptr = (int*) map;
	for (i = 0; i < num_values; i++) {
		zipf_rv = zipf(alpha, n);
		if (norm == 1) {
			int pr = (int) ceil(log10(n));
			double zipf_float = ((double) zipf_rv / n);
			iptr[0] = pr;
			fptr = (float*) (iptr + 1);
			fptr[0] = zipf_float;
			iptr = (int*) (fptr + 1);
		} else {
			int zipf_int = (int) zipf_rv;
			float r = ((double) rand() / (RAND_MAX));
			iptr[0] = zipf_int;
			fptr = (float*) iptr + 1;
			fptr[0] = r;
			iptr = (int*) fptr + 1;
			//fprintf(fp, "%d\t%.10F\n", zipf_int, r);
		}
	}

	// Output "done" message and close the output file
	printf("-------------------------------------------------------- \n");
	printf("-  Done! \n");
	printf("-------------------------------------------------------- \n");

	if (munmap(map, FILESIZE) == -1) {
		perror("Error un-mmapping the file");
		/* Decide here whether to close(fd) and exit() or not. Depends... */
	}

	close(fd);
	return -1;
}

//===========================================================================
//=  Function to generate Zipf (power law) distributed random variables     =
//=    - Input: alpha and N                                                 =
//=    - Output: Returns with Zipf distributed random variable              =
//===========================================================================
int zipf(double alpha, int n) {
	static int first = TRUE;      // Static first time flag
	static double c = 0;          // Normalization constant
	double z;                     // Uniform random number (0 < z < 1)
	double sum_prob;              // Sum of probabilities
	double zipf_value;            // Computed exponential value to be returned
	int i;                     // Loop counter

	// Compute normalization constant on first call only
	if (first == TRUE) {
		for (i = 1; i <= n; i++)
			c = c + (1.0 / pow((double) i, alpha));
		c = 1.0 / c;
		first = FALSE;
	}

	// Pull a uniform random number (0 < z < 1)
	do {
		z = rand_val(0);
	} while ((z == 0) || (z == 1));

	// Map z to the value
	sum_prob = 0;
	for (i = 1; i <= n; i++) {
		sum_prob = sum_prob + c / pow((double) i, alpha);
		if (sum_prob >= z) {
			zipf_value = i;
			break;
		}
	}

	// Assert that zipf_value is between 1 and N
	assert((zipf_value >= 1) && (zipf_value <= n));

	return (zipf_value);
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=========================================================================
double rand_val(int seed) {
	const long a = 16807;  // Multiplier
	const long m = 2147483647;  // Modulus
	const long q = 127773;  // m div a
	const long r = 2836;  // m mod a
	static long x;               // Random int value
	long x_div_q;         // x divided by q
	long x_mod_q;         // x modulo q
	long x_new;           // New x value

	// Set the seed if argument is non-zero and then return zero
	if (seed > 0) {
		x = seed;
		return (0.0);
	}

	// RNG using integer arithmetic
	x_div_q = x / q;
	x_mod_q = x % q;
	x_new = (a * x_mod_q) - (r * x_div_q);
	if (x_new > 0)
		x = x_new;
	else
		x = x_new + m;

	// Return a random value between 0.0 and 1.0
	return ((double) x / m);
}
